!--------------------------------------------------------------------------------------------------!
!   CP2K: A general program to perform molecular dynamics simulations                              !
!   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
!                                                                                                  !
!   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
!--------------------------------------------------------------------------------------------------!

! **************************************************************************************************
!> \brief contains a functional that calculates the energy and its derivatives
!>      for the geometry optimizer
!> \par History
!>      none
! **************************************************************************************************
MODULE gopt_f_methods

   USE atomic_kind_list_types, ONLY: atomic_kind_list_type
   USE atomic_kind_types, ONLY: atomic_kind_type, &
                                get_atomic_kind_set
   USE cell_methods, ONLY: cell_create, &
                           init_cell
   USE cell_types, ONLY: cell_copy, &
                         cell_release, &
                         cell_type, &
                         real_to_scaled, &
                         scaled_to_real
   USE cp_log_handling, ONLY: cp_to_string
   USE cp_subsys_types, ONLY: cp_subsys_get, &
                              cp_subsys_set, &
                              cp_subsys_type, &
                              pack_subsys_particles
   USE cp_units, ONLY: cp_unit_from_cp2k
   USE dimer_types, ONLY: dimer_env_type
   USE dimer_utils, ONLY: update_dimer_vec
   USE distribution_1d_types, ONLY: distribution_1d_type
   USE force_env_types, ONLY: force_env_get, &
                              force_env_get_natom, &
                              force_env_get_nparticle, &
                              force_env_type, &
                              use_qmmm, &
                              use_qmmmx
   USE gopt_f_types, ONLY: gopt_f_type
   USE gopt_param_types, ONLY: gopt_param_type
   USE input_constants, ONLY: default_cell_direct_id, &
                              default_cell_geo_opt_id, &
                              default_cell_md_id, &
                              default_cell_method_id, &
                              default_minimization_method_id, &
                              default_shellcore_method_id, &
                              default_ts_method_id, &
                              fix_none
   USE input_cp2k_restarts, ONLY: write_restart
   USE input_section_types, ONLY: section_vals_type, &
                                  section_vals_val_get
   USE kinds, ONLY: default_string_length, &
                    dp, &
                    int_8
   USE machine, ONLY: m_flush
   USE md_energies, ONLY: sample_memory
   USE message_passing, ONLY: mp_para_env_type
   USE motion_utils, ONLY: write_simulation_cell, &
                           write_stress_tensor, &
                           write_trajectory
   USE particle_list_types, ONLY: particle_list_type
   USE particle_methods, ONLY: write_structure_data
   USE particle_types, ONLY: particle_type
   USE qmmm_util, ONLY: apply_qmmm_translate
   USE qmmmx_util, ONLY: apply_qmmmx_translate
   USE virial_methods, ONLY: virial_evaluate
   USE virial_types, ONLY: virial_type
#include "../base/base_uses.f90"

   IMPLICIT NONE
   PRIVATE

   #:include "gopt_f77_methods.fypp"

   LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'gopt_f_methods'

   PUBLIC :: gopt_f_create_x0, &
             print_geo_opt_header, print_geo_opt_nc, &
             gopt_f_io_init, gopt_f_io, gopt_f_io_finalize, gopt_f_ii, &
             apply_cell_change

CONTAINS

! **************************************************************************************************
!> \brief returns the value of the parameters for the actual configuration
!> \param gopt_env the geometry optimization environment you want the info about
!>      x0: the parameter vector (is allocated by this routine)
!> \param x0 ...
!> \par History
!>      - Cell optimization revised (06.11.2012,MK)
! **************************************************************************************************
   SUBROUTINE gopt_f_create_x0(gopt_env, x0)

      TYPE(gopt_f_type), POINTER                         :: gopt_env
      REAL(KIND=dp), DIMENSION(:), POINTER               :: x0

      INTEGER                                            :: i, idg, j, nparticle
      TYPE(cell_type), POINTER                           :: cell
      TYPE(cp_subsys_type), POINTER                      :: subsys

      NULLIFY (cell)
      NULLIFY (subsys)

      SELECT CASE (gopt_env%type_id)
      CASE (default_minimization_method_id, default_ts_method_id)
         CALL force_env_get(gopt_env%force_env, subsys=subsys)
         ! before starting we handle the case of translating coordinates (QM/MM)
         IF (gopt_env%force_env%in_use == use_qmmm) &
            CALL apply_qmmm_translate(gopt_env%force_env%qmmm_env)
         IF (gopt_env%force_env%in_use == use_qmmmx) &
            CALL apply_qmmmx_translate(gopt_env%force_env%qmmmx_env)
         nparticle = force_env_get_nparticle(gopt_env%force_env)
         ALLOCATE (x0(3*nparticle))
         CALL pack_subsys_particles(subsys=subsys, r=x0)
      CASE (default_cell_method_id)
         SELECT CASE (gopt_env%cell_method_id)
         CASE (default_cell_direct_id)
            CALL force_env_get(gopt_env%force_env, subsys=subsys, cell=cell)
            ! Store reference cell
            gopt_env%h_ref = cell%hmat
            ! before starting we handle the case of translating coordinates (QM/MM)
            IF (gopt_env%force_env%in_use == use_qmmm) &
               CALL apply_qmmm_translate(gopt_env%force_env%qmmm_env)
            IF (gopt_env%force_env%in_use == use_qmmmx) &
               CALL apply_qmmmx_translate(gopt_env%force_env%qmmmx_env)
            nparticle = force_env_get_nparticle(gopt_env%force_env)
            ALLOCATE (x0(3*nparticle + 6))
            CALL pack_subsys_particles(subsys=subsys, r=x0)
            idg = 3*nparticle
            DO i = 1, 3
               DO j = 1, i
                  idg = idg + 1
                  x0(idg) = cell%hmat(j, i)
               END DO
            END DO
         CASE (default_cell_geo_opt_id, default_cell_md_id)
            CALL force_env_get(gopt_env%force_env, cell=cell)
            ALLOCATE (x0(6))
            idg = 0
            DO i = 1, 3
               DO j = 1, i
                  idg = idg + 1
                  x0(idg) = cell%hmat(j, i)
               END DO
            END DO
         CASE DEFAULT
            CPABORT("")
         END SELECT
      CASE DEFAULT
         CPABORT("")
      END SELECT

   END SUBROUTINE gopt_f_create_x0

! **************************************************************************************************
!> \brief Prints iteration step of the optimization procedure on screen
!> \param its ...
!> \param output_unit ...
!> \author Teodoro Laino [tlaino] - University of Zurich - 03.2008
! **************************************************************************************************
   SUBROUTINE gopt_f_ii(its, output_unit)

      INTEGER, INTENT(IN)                                :: its, output_unit

      IF (output_unit > 0) THEN
         WRITE (UNIT=output_unit, FMT="(/,T2,26('-'))")
         WRITE (UNIT=output_unit, FMT="(T2,A,I6)") "OPTIMIZATION STEP: ", its
         WRITE (UNIT=output_unit, FMT="(T2,26('-'))")
         CALL m_flush(output_unit)
      END IF

   END SUBROUTINE gopt_f_ii

! **************************************************************************************************
!> \brief Handles the Output during an optimization run
!> \param gopt_env ...
!> \param output_unit ...
!> \param opt_energy ...
!> \param wildcard ...
!> \param its ...
!> \param used_time ...
!> \author Teodoro Laino [tlaino] - University of Zurich - 03.2008
! **************************************************************************************************
   SUBROUTINE gopt_f_io_init(gopt_env, output_unit, opt_energy, wildcard, its, used_time)
      TYPE(gopt_f_type), POINTER                         :: gopt_env
      INTEGER, INTENT(IN)                                :: output_unit
      REAL(KIND=dp)                                      :: opt_energy
      CHARACTER(LEN=5)                                   :: wildcard
      INTEGER, INTENT(IN)                                :: its
      REAL(KIND=dp)                                      :: used_time

      TYPE(mp_para_env_type), POINTER                    :: para_env
      REAL(KIND=dp)                                      :: pres_int
      INTEGER(KIND=int_8)                                :: max_memory
      LOGICAL                                            :: print_memory

      NULLIFY (para_env)
      CALL section_vals_val_get(gopt_env%motion_section, "PRINT%MEMORY_INFO", l_val=print_memory)
      max_memory = 0
      IF (print_memory) THEN
         CALL force_env_get(gopt_env%force_env, para_env=para_env)
         max_memory = sample_memory(para_env)
      END IF

      SELECT CASE (gopt_env%type_id)
      CASE (default_ts_method_id, default_minimization_method_id)
         ! Geometry Optimization (Minimization and Transition State Search)
         IF (.NOT. gopt_env%dimer_rotation) THEN
            CALL write_cycle_infos(output_unit, it=its, etot=opt_energy, wildcard=wildcard, &
                                   used_time=used_time, max_memory=max_memory)
         ELSE
            CALL write_rot_cycle_infos(output_unit, it=its, etot=opt_energy, dimer_env=gopt_env%dimer_env, &
                                       wildcard=wildcard, used_time=used_time, max_memory=max_memory)
         END IF
      CASE (default_cell_method_id)
         ! Cell Optimization
         pres_int = gopt_env%cell_env%pres_int
         CALL write_cycle_infos(output_unit, it=its, etot=opt_energy, pres_int=pres_int, &
                                wildcard=wildcard, used_time=used_time, max_memory=max_memory)
      CASE (default_shellcore_method_id)
         CALL write_cycle_infos(output_unit, it=its, etot=opt_energy, wildcard=wildcard, &
                                used_time=used_time, max_memory=max_memory)
      END SELECT

   END SUBROUTINE gopt_f_io_init

! **************************************************************************************************
!> \brief Handles the Output during an optimization run
!> \param gopt_env ...
!> \param force_env ...
!> \param root_section ...
!> \param its ...
!> \param opt_energy ...
!> \param output_unit ...
!> \param eold ...
!> \param emin ...
!> \param wildcard ...
!> \param gopt_param ...
!> \param ndf ...
!> \param dx ...
!> \param xi ...
!> \param conv ...
!> \param pred ...
!> \param rat ...
!> \param step ...
!> \param rad ...
!> \param used_time ...
!> \author Teodoro Laino [tlaino] - University of Zurich - 03.2008
! **************************************************************************************************
   SUBROUTINE gopt_f_io(gopt_env, force_env, root_section, its, opt_energy, &
                        output_unit, eold, emin, wildcard, gopt_param, ndf, dx, xi, conv, pred, rat, &
                        step, rad, used_time)
      TYPE(gopt_f_type), POINTER                         :: gopt_env
      TYPE(force_env_type), POINTER                      :: force_env
      TYPE(section_vals_type), POINTER                   :: root_section
      INTEGER, INTENT(IN)                                :: its
      REAL(KIND=dp), INTENT(IN)                          :: opt_energy
      INTEGER, INTENT(IN)                                :: output_unit
      REAL(KIND=dp)                                      :: eold, emin
      CHARACTER(LEN=5)                                   :: wildcard
      TYPE(gopt_param_type), POINTER                     :: gopt_param
      INTEGER, INTENT(IN), OPTIONAL                      :: ndf
      REAL(KIND=dp), DIMENSION(:), INTENT(IN), OPTIONAL  :: dx
      REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER     :: xi
      LOGICAL, OPTIONAL                                  :: conv
      REAL(KIND=dp), INTENT(IN), OPTIONAL                :: pred, rat, step, rad
      REAL(KIND=dp)                                      :: used_time

      INTEGER(KIND=int_8)                                :: max_memory
      LOGICAL                                            :: print_memory
      REAL(KIND=dp)                                      :: pres_diff, pres_diff_constr, pres_int, &
                                                            pres_tol
      TYPE(mp_para_env_type), POINTER                    :: para_env

      NULLIFY (para_env)
      CALL section_vals_val_get(gopt_env%motion_section, "PRINT%MEMORY_INFO", l_val=print_memory)
      max_memory = 0
      IF (print_memory) THEN
         CALL force_env_get(force_env, para_env=para_env)
         max_memory = sample_memory(para_env)
      END IF

      SELECT CASE (gopt_env%type_id)
      CASE (default_ts_method_id, default_minimization_method_id)
         ! Geometry Optimization (Minimization and Transition State Search)
         IF (.NOT. gopt_env%dimer_rotation) THEN
            CALL geo_opt_io(force_env=force_env, root_section=root_section, &
                            motion_section=gopt_env%motion_section, its=its, opt_energy=opt_energy)
            CALL write_cycle_infos(output_unit, its, etot=opt_energy, ediff=opt_energy - eold, &
                                   pred=pred, rat=rat, step=step, rad=rad, emin=emin, &
                                   wildcard=wildcard, used_time=used_time, max_memory=max_memory)
            ! Possibly check convergence
            IF (PRESENT(conv)) THEN
               CPASSERT(PRESENT(ndf))
               CPASSERT(PRESENT(dx))
               CPASSERT(PRESENT(xi))
               CALL check_converg(ndf, dx, xi, output_unit, conv, gopt_param, max_memory)
            END IF
         ELSE
            CALL update_dimer_vec(gopt_env%dimer_env, gopt_env%motion_section)
            CALL write_restart(force_env=force_env, root_section=root_section)
            CALL write_rot_cycle_infos(output_unit, its, opt_energy, opt_energy - eold, emin, gopt_env%dimer_env, &
                                       wildcard=wildcard, used_time=used_time, max_memory=max_memory)
            ! Possibly check convergence
            IF (PRESENT(conv)) THEN
               CPASSERT(ASSOCIATED(gopt_env%dimer_env))
               CALL check_rot_conv(gopt_env%dimer_env, output_unit, conv)
            END IF
         END IF
      CASE (default_cell_method_id)
         ! Cell Optimization
         pres_diff = gopt_env%cell_env%pres_int - gopt_env%cell_env%pres_ext
         pres_int = gopt_env%cell_env%pres_int
         pres_tol = gopt_env%cell_env%pres_tol
         CALL geo_opt_io(force_env=force_env, root_section=root_section, &
                         motion_section=gopt_env%motion_section, its=its, opt_energy=opt_energy)
         CALL write_cycle_infos(output_unit, its, etot=opt_energy, ediff=opt_energy - eold, &
                                pred=pred, rat=rat, step=step, rad=rad, emin=emin, pres_int=pres_int, &
                                wildcard=wildcard, used_time=used_time, max_memory=max_memory)
         ! Possibly check convergence
         IF (PRESENT(conv)) THEN
            CPASSERT(PRESENT(ndf))
            CPASSERT(PRESENT(dx))
            CPASSERT(PRESENT(xi))
            IF (gopt_env%cell_env%constraint_id == fix_none) THEN
               CALL check_converg(ndf, dx, xi, output_unit, conv, gopt_param, max_memory, pres_diff, pres_tol)
            ELSE
               pres_diff_constr = gopt_env%cell_env%pres_constr - gopt_env%cell_env%pres_ext
               CALL check_converg(ndf, dx, xi, output_unit, conv, gopt_param, max_memory, pres_diff, &
                                  pres_tol, pres_diff_constr)
            END IF
         END IF
      CASE (default_shellcore_method_id)
         CALL write_cycle_infos(output_unit, its, etot=opt_energy, ediff=opt_energy - eold, &
                                pred=pred, rat=rat, step=step, rad=rad, emin=emin, wildcard=wildcard, &
                                used_time=used_time, max_memory=max_memory)
         ! Possibly check convergence
         IF (PRESENT(conv)) THEN
            CPASSERT(PRESENT(ndf))
            CPASSERT(PRESENT(dx))
            CPASSERT(PRESENT(xi))
            CALL check_converg(ndf, dx, xi, output_unit, conv, gopt_param, max_memory)
         END IF
      END SELECT

   END SUBROUTINE gopt_f_io

! **************************************************************************************************
!> \brief Handles the Output at the end of an optimization run
!> \param gopt_env ...
!> \param force_env ...
!> \param x0 ...
!> \param conv ...
!> \param its ...
!> \param root_section ...
!> \param para_env ...
!> \param master ...
!> \param output_unit ...
!> \author Teodoro Laino [tlaino] - University of Zurich - 03.2008
! **************************************************************************************************
   RECURSIVE SUBROUTINE gopt_f_io_finalize(gopt_env, force_env, x0, conv, its, root_section, &
                                           para_env, master, output_unit)
      TYPE(gopt_f_type), POINTER                         :: gopt_env
      TYPE(force_env_type), POINTER                      :: force_env
      REAL(KIND=dp), DIMENSION(:), POINTER               :: x0
      LOGICAL                                            :: conv
      INTEGER                                            :: its
      TYPE(section_vals_type), POINTER                   :: root_section
      TYPE(mp_para_env_type), POINTER                    :: para_env
      INTEGER, INTENT(IN)                                :: master, output_unit

      IF (gopt_env%eval_opt_geo) THEN
         IF (.NOT. gopt_env%dimer_rotation) THEN
            CALL write_final_info(output_unit, conv, its, gopt_env, x0, master, &
                                  para_env, force_env, gopt_env%motion_section, root_section)
         ELSE
            CALL update_dimer_vec(gopt_env%dimer_env, gopt_env%motion_section)
            CALL write_restart(force_env=force_env, root_section=root_section)
         END IF
      END IF

   END SUBROUTINE gopt_f_io_finalize

! **************************************************************************************************
!> \brief ...
!> \param output_unit ...
!> \param it ...
!> \param etot ...
!> \param ediff ...
!> \param pred ...
!> \param rat ...
!> \param step ...
!> \param rad ...
!> \param emin ...
!> \param pres_int ...
!> \param wildcard ...
!> \param used_time ...
! **************************************************************************************************
   SUBROUTINE write_cycle_infos(output_unit, it, etot, ediff, pred, rat, step, rad, emin, &
                                pres_int, wildcard, used_time, max_memory)

      INTEGER, INTENT(IN)                                :: output_unit, it
      REAL(KIND=dp), INTENT(IN)                          :: etot
      REAL(KIND=dp), INTENT(IN), OPTIONAL                :: ediff, pred, rat, step, rad, emin, &
                                                            pres_int
      CHARACTER(LEN=5), INTENT(IN)                       :: wildcard
      REAL(KIND=dp), INTENT(IN)                          :: used_time
      INTEGER(KIND=int_8), INTENT(IN)                    :: max_memory

      CHARACTER(LEN=5)                                   :: tag
      REAL(KIND=dp)                                      :: tmp_r1

      IF (output_unit > 0) THEN
         tag = "OPT| "
         WRITE (UNIT=output_unit, FMT="(/,T2,A)") tag//REPEAT("*", 74)
         WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,I25)") &
            tag//"Step number", it
         WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,A25)") &
            tag//"Optimization method", wildcard
         WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,F25.10)") &
            tag//"Total energy [hartree]", etot
         IF (PRESENT(pres_int)) THEN
            tmp_r1 = cp_unit_from_cp2k(pres_int, "bar")
            WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,F25.10)") &
               tag//"Internal pressure [bar]", tmp_r1
         END IF
         IF (PRESENT(ediff)) THEN
            WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,F25.10)") &
               tag//"Effective energy change [hartree]", ediff
         END IF
         IF (PRESENT(pred)) THEN
            WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,F25.10)") &
               tag//"Predicted energy change [hartree]", pred
         END IF
         IF (PRESENT(rat)) THEN
            WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,F25.10)") &
               tag//"Scaling factor", rat
         END IF
         IF (PRESENT(step)) THEN
            WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,F25.10)") &
               tag//"Step size", step
         END IF
         IF (PRESENT(rad)) THEN
            WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,F25.10)") &
               tag//"Trust radius", rad
         END IF
         IF (PRESENT(emin)) THEN
            IF (etot < emin) THEN
               WRITE (UNIT=output_unit, FMT="(T2,A,T77,A4)") &
                  tag//"Decrease in energy", " YES"
            ELSE
               WRITE (UNIT=output_unit, FMT="(T2,A,T77,A4)") &
                  tag//"Decrease in energy", "  NO"
            END IF
         END IF
         WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,F25.3)") &
            tag//"Used time [s]", used_time
         IF (it == 0) THEN
            WRITE (UNIT=output_unit, FMT="(T2,A)") tag//REPEAT("*", 74)
            IF (max_memory /= 0) THEN
               WRITE (UNIT=output_unit, FMT="(T2,A,T60,1X,I20)") &
                  tag//"Estimated peak process memory [MiB]", &
                  (max_memory + (1024*1024) - 1)/(1024*1024)
            END IF
         END IF
      END IF

   END SUBROUTINE write_cycle_infos

! **************************************************************************************************
!> \brief ...
!> \param output_unit ...
!> \param it ...
!> \param etot ...
!> \param ediff ...
!> \param emin ...
!> \param dimer_env ...
!> \param used_time ...
!> \param wildcard ...
!> \date  01.2008
!> \author Luca Bellucci and Teodoro Laino - created [tlaino]
! **************************************************************************************************
   SUBROUTINE write_rot_cycle_infos(output_unit, it, etot, ediff, emin, dimer_env, used_time, &
                                    wildcard, max_memory)

      INTEGER, INTENT(IN)                                :: output_unit, it
      REAL(KIND=dp), INTENT(IN)                          :: etot
      REAL(KIND=dp), INTENT(IN), OPTIONAL                :: ediff, emin
      TYPE(dimer_env_type), POINTER                      :: dimer_env
      REAL(KIND=dp), INTENT(IN)                          :: used_time
      CHARACTER(LEN=5), INTENT(IN)                       :: wildcard
      INTEGER(KIND=int_8), INTENT(IN)                    :: max_memory

      CHARACTER(LEN=5)                                   :: tag

      IF (output_unit > 0) THEN
         tag = "OPT| "
         WRITE (UNIT=output_unit, FMT="(/,T2,A)") tag//REPEAT("*", 74)
         WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,I25)") &
            tag//"Rotational step number", it
         WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,A25)") &
            tag//"Optimization method", wildcard
         WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,F25.10)") &
            tag//"Local curvature", dimer_env%rot%curvature, &
            tag//"Total rotational force", etot
         IF (PRESENT(ediff)) THEN
            WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,F25.10)") &
               tag//"Rotational force change", ediff
         END IF
         IF (PRESENT(emin)) THEN
            IF (etot < emin) THEN
               WRITE (UNIT=output_unit, FMT="(T2,A,T77,A4)") &
                  tag//"Decrease in rotational force", " YES"
            ELSE
               WRITE (UNIT=output_unit, FMT="(T2,A,T77,A4)") &
                  tag//"Decrease in rotational force", "  NO"
            END IF
         END IF
         WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,F25.3)") &
            tag//"Used time [s]", used_time
         IF (it == 0) THEN
            WRITE (UNIT=output_unit, FMT="(T2,A)") tag//REPEAT("*", 74)
            IF (max_memory /= 0) THEN
               WRITE (UNIT=output_unit, FMT="(T2,A,T60,1X,I20)") &
                  tag//"Estimated peak process memory [MiB]", &
                  (max_memory + (1024*1024) - 1)/(1024*1024)
            END IF
         END IF
      END IF

   END SUBROUTINE write_rot_cycle_infos

! **************************************************************************************************
!> \brief ...
!> \param ndf ...
!> \param dr ...
!> \param g ...
!> \param output_unit ...
!> \param conv ...
!> \param gopt_param ...
!> \param max_memory ...
!> \param pres_diff ...
!> \param pres_tol ...
!> \param pres_diff_constr ...
! **************************************************************************************************
   SUBROUTINE check_converg(ndf, dr, g, output_unit, conv, gopt_param, max_memory, pres_diff, &
                            pres_tol, pres_diff_constr)

      INTEGER, INTENT(IN)                                :: ndf
      REAL(KIND=dp), INTENT(IN)                          :: dr(ndf), g(ndf)
      INTEGER, INTENT(IN)                                :: output_unit
      LOGICAL, INTENT(OUT)                               :: conv
      TYPE(gopt_param_type), POINTER                     :: gopt_param
      INTEGER(KIND=int_8), INTENT(IN)                    :: max_memory
      REAL(KIND=dp), INTENT(IN), OPTIONAL                :: pres_diff, pres_tol, pres_diff_constr

      CHARACTER(LEN=5)                                   :: tag
      INTEGER                                            :: indf
      LOGICAL                                            :: conv_dx, conv_g, conv_p, conv_rdx, &
                                                            conv_rg
      REAL(KIND=dp)                                      :: dumm, dxcon, gcon, maxdum(4), rmsgcon, &
                                                            rmsxcon, tmp_r1

      dxcon = gopt_param%max_dr
      gcon = gopt_param%max_force
      rmsgcon = gopt_param%rms_force
      rmsxcon = gopt_param%rms_dr

      conv = .FALSE.
      conv_dx = .TRUE.
      conv_rdx = .TRUE.
      conv_g = .TRUE.
      conv_rg = .TRUE.
      conv_p = .TRUE.

      dumm = 0.0_dp
      DO indf = 1, ndf
         IF (indf == 1) maxdum(1) = ABS(dr(indf))
         dumm = dumm + dr(indf)**2
         IF (ABS(dr(indf)) > dxcon) conv_dx = .FALSE.
         IF (ABS(dr(indf)) > maxdum(1)) maxdum(1) = ABS(dr(indf))
      END DO
      ! SQRT(dumm/ndf) > rmsxcon
      IF (dumm > (rmsxcon*rmsxcon*ndf)) conv_rdx = .FALSE.
      maxdum(2) = SQRT(dumm/ndf)

      dumm = 0.0_dp
      DO indf = 1, ndf
         IF (indf == 1) maxdum(3) = ABS(g(indf))
         dumm = dumm + g(indf)**2
         IF (ABS(g(indf)) > gcon) conv_g = .FALSE.
         IF (ABS(g(indf)) > maxdum(3)) maxdum(3) = ABS(g(indf))
      END DO
      ! SQRT(dumm/ndf) > rmsgcon
      IF (dumm > (rmsgcon*rmsgcon*ndf)) conv_rg = .FALSE.
      maxdum(4) = SQRT(dumm/ndf)

      IF (PRESENT(pres_diff_constr) .AND. PRESENT(pres_tol)) THEN
         conv_p = ABS(pres_diff_constr) < ABS(pres_tol)
      ELSEIF (PRESENT(pres_diff) .AND. PRESENT(pres_tol)) THEN
         conv_p = ABS(pres_diff) < ABS(pres_tol)
      END IF

      IF (output_unit > 0) THEN

         tag = "OPT| "

         WRITE (UNIT=output_unit, FMT="(T2,A)") TRIM(tag)
         WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,F25.10)") &
            tag//"Maximum step size", maxdum(1), &
            tag//"Convergence limit for maximum step size", dxcon
         IF (conv_dx) THEN
            WRITE (UNIT=output_unit, FMT="(T2,A,T77,A4)") &
               tag//"Maximum step size is converged", " YES"
         ELSE
            WRITE (UNIT=output_unit, FMT="(T2,A,T77,A4)") &
               tag//"Maximum step size is converged", "  NO"
         END IF

         WRITE (UNIT=output_unit, FMT="(T2,A)") TRIM(tag)
         WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,F25.10)") &
            tag//"RMS step size", maxdum(2), &
            tag//"Convergence limit for RMS step size", rmsxcon
         IF (conv_rdx) THEN
            WRITE (UNIT=output_unit, FMT="(T2,A,T77,A4)") &
               tag//"RMS step size is converged", " YES"
         ELSE
            WRITE (UNIT=output_unit, FMT="(T2,A,T77,A4)") &
               tag//"RMS step size is converged", "  NO"
         END IF

         WRITE (UNIT=output_unit, FMT="(T2,A)") TRIM(tag)
         WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,F25.10)") &
            tag//"Maximum gradient", maxdum(3), &
            tag//"Convergence limit for maximum gradient", gcon
         IF (conv_g) THEN
            WRITE (UNIT=output_unit, FMT="(T2,A,T77,A4)") &
               tag//"Maximum gradient is converged", " YES"
         ELSE
            WRITE (UNIT=output_unit, FMT="(T2,A,T77,A4)") &
               tag//"Maximum gradient is converged", "  NO"
         END IF

         WRITE (UNIT=output_unit, FMT="(T2,A)") TRIM(tag)
         WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,F25.10)") &
            tag//"RMS gradient", maxdum(4), &
            tag//"Convergence limit for RMS gradient", rmsgcon
         IF (conv_rg) THEN
            WRITE (UNIT=output_unit, FMT="(T2,A,T77,A4)") &
               tag//"RMS gradient is converged", " YES"
         ELSE
            WRITE (UNIT=output_unit, FMT="(T2,A,T77,A4)") &
               tag//"RMS gradient is converged", "  NO"
         END IF

         IF (PRESENT(pres_diff) .AND. PRESENT(pres_tol)) THEN
            tmp_r1 = cp_unit_from_cp2k(pres_diff, "bar")
            WRITE (UNIT=output_unit, FMT="(T2,A)") TRIM(tag)
            IF (PRESENT(pres_diff_constr)) THEN
               WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,F25.10)") &
                  tag//"Pressure deviation without constraint [bar]", tmp_r1
               tmp_r1 = cp_unit_from_cp2k(pres_diff_constr, "bar")
               WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,F25.10)") &
                  tag//"Pressure deviation with constraint [bar]", tmp_r1
            ELSE
               WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,F25.10)") &
                  tag//"Pressure deviation [bar]", tmp_r1
            END IF
            tmp_r1 = cp_unit_from_cp2k(pres_tol, "bar")
            WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,F25.10)") &
               tag//"Pressure tolerance [bar]", tmp_r1
            IF (conv_p) THEN
               WRITE (UNIT=output_unit, FMT="(T2,A,T77,A4)") &
                  tag//"Pressure is converged", " YES"
            ELSE
               WRITE (UNIT=output_unit, FMT="(T2,A,T77,A4)") &
                  tag//"Pressure is converged", "  NO"
            END IF
         END IF

         WRITE (UNIT=output_unit, FMT="(T2,A)") tag//REPEAT("*", 74)

         IF (max_memory /= 0) THEN
            WRITE (UNIT=output_unit, FMT="(T2,A,T60,1X,I20)") &
               tag//"Estimated peak process memory after this step [MiB]", &
               (max_memory + (1024*1024) - 1)/(1024*1024)
         END IF

      END IF

      IF (conv_dx .AND. conv_rdx .AND. conv_g .AND. conv_rg .AND. conv_p) conv = .TRUE.

      IF ((conv) .AND. (output_unit > 0)) THEN
         WRITE (UNIT=output_unit, FMT="(/,T2,A)") REPEAT("*", 79)
         WRITE (UNIT=output_unit, FMT="(T2,A,T25,A,T78,A)") &
            "***", "GEOMETRY OPTIMIZATION COMPLETED", "***"
         WRITE (UNIT=output_unit, FMT="(T2,A)") REPEAT("*", 79)
      END IF

   END SUBROUTINE check_converg

! **************************************************************************************************
!> \brief ...
!> \param dimer_env ...
!> \param output_unit ...
!> \param conv ...
!> \date  01.2008
!> \author Luca Bellucci and Teodoro Laino - created [tlaino]
! **************************************************************************************************
   SUBROUTINE check_rot_conv(dimer_env, output_unit, conv)

      TYPE(dimer_env_type), POINTER                      :: dimer_env
      INTEGER, INTENT(IN)                                :: output_unit
      LOGICAL, INTENT(OUT)                               :: conv

      CHARACTER(LEN=5)                                   :: tag

      conv = (ABS(dimer_env%rot%angle2) < dimer_env%rot%angle_tol)

      IF (output_unit > 0) THEN
         tag = "OPT| "
         WRITE (UNIT=output_unit, FMT="(T2,A)") TRIM(tag)
         WRITE (UNIT=output_unit, FMT="(T2,A,T55,1X,F25.10)") &
            tag//"Predicted angle step size", dimer_env%rot%angle1, &
            tag//"Effective angle step size", dimer_env%rot%angle2, &
            tag//"Convergence limit for angle step size", dimer_env%rot%angle_tol
         IF (conv) THEN
            WRITE (UNIT=output_unit, FMT="(T2,A,T77,A4)") &
               tag//"Angle step size is converged", " YES"
         ELSE
            WRITE (UNIT=output_unit, FMT="(T2,A,T77,A4)") &
               tag//"Angle step size is converged", "  NO"
         END IF
         WRITE (UNIT=output_unit, FMT="(T2,A)") tag//REPEAT("*", 74)
      END IF

      IF ((conv) .AND. (output_unit > 0)) THEN
         WRITE (UNIT=output_unit, FMT="(/,T2,A)") REPEAT("*", 79)
         WRITE (UNIT=output_unit, FMT="(T2,A,T25,A,T78,A)") &
            "***", "ROTATION OPTIMIZATION COMPLETED", "***"
         WRITE (UNIT=output_unit, FMT="(T2,A)") REPEAT("*", 79)
      END IF

   END SUBROUTINE check_rot_conv

! **************************************************************************************************
!> \brief ...
!> \param output_unit ...
!> \param conv ...
!> \param it ...
!> \param gopt_env ...
!> \param x0 ...
!> \param master ...
!> \param para_env ...
!> \param force_env ...
!> \param motion_section ...
!> \param root_section ...
!> \date  11.2007
!> \author Teodoro Laino [tlaino] - University of Zurich
! **************************************************************************************************
   RECURSIVE SUBROUTINE write_final_info(output_unit, conv, it, gopt_env, x0, master, para_env, force_env, &
                                         motion_section, root_section)
      INTEGER, INTENT(IN)                                :: output_unit
      LOGICAL, INTENT(IN)                                :: conv
      INTEGER, INTENT(INOUT)                             :: it
      TYPE(gopt_f_type), POINTER                         :: gopt_env
      REAL(KIND=dp), DIMENSION(:), POINTER               :: x0
      INTEGER, INTENT(IN)                                :: master
      TYPE(mp_para_env_type), POINTER                    :: para_env
      TYPE(force_env_type), POINTER                      :: force_env
      TYPE(section_vals_type), POINTER                   :: motion_section, root_section

      REAL(KIND=dp)                                      :: etot
      TYPE(cell_type), POINTER                           :: cell
      TYPE(cp_subsys_type), POINTER                      :: subsys
      TYPE(particle_list_type), POINTER                  :: particles
      TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set

      CALL force_env_get(force_env, cell=cell, subsys=subsys)
      CALL cp_subsys_get(subsys=subsys, particles=particles)
      particle_set => particles%els
      IF (conv) THEN
         it = it + 1
         CALL write_structure_data(particle_set, cell, motion_section)
         CALL write_restart(force_env=force_env, root_section=root_section)

         IF (output_unit > 0) &
            WRITE (UNIT=output_unit, FMT="(/,T20,' Reevaluating energy at the minimum')")

         CALL cp_eval_at(gopt_env, x0, f=etot, master=master, final_evaluation=.TRUE., &
                         para_env=para_env)
         CALL write_geo_traj(force_env, root_section, it, etot)
      END IF

   END SUBROUTINE write_final_info

! **************************************************************************************************
!> \brief  Specific driver for dumping trajectory during a GEO_OPT
!> \param force_env ...
!> \param root_section ...
!> \param it ...
!> \param etot ...
!> \date   11.2007
!> \par    History
!>         09.2010: Output of core and shell positions and forces (MK)
!> \author Teodoro Laino [tlaino] - University of Zurich
! **************************************************************************************************
   SUBROUTINE write_geo_traj(force_env, root_section, it, etot)

      TYPE(force_env_type), POINTER                      :: force_env
      TYPE(section_vals_type), POINTER                   :: root_section
      INTEGER, INTENT(IN)                                :: it
      REAL(KIND=dp), INTENT(IN)                          :: etot

      LOGICAL                                            :: shell_adiabatic, shell_present
      TYPE(atomic_kind_list_type), POINTER               :: atomic_kinds
      TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
      TYPE(cp_subsys_type), POINTER                      :: subsys
      TYPE(particle_list_type), POINTER                  :: core_particles, shell_particles

      NULLIFY (atomic_kinds)
      NULLIFY (atomic_kind_set)
      NULLIFY (core_particles)
      NULLIFY (shell_particles)
      NULLIFY (subsys)

      CALL write_trajectory(force_env, root_section, it, 0.0_dp, 0.0_dp, etot)
      ! Print Force
      CALL write_trajectory(force_env, root_section, it, 0.0_dp, 0.0_dp, etot, "FORCES", middle_name="frc")
      CALL force_env_get(force_env, subsys=subsys)
      CALL cp_subsys_get(subsys, atomic_kinds=atomic_kinds)
      atomic_kind_set => atomic_kinds%els
      CALL get_atomic_kind_set(atomic_kind_set, &
                               shell_present=shell_present, &
                               shell_adiabatic=shell_adiabatic)
      IF (shell_present) THEN
         CALL cp_subsys_get(subsys, &
                            core_particles=core_particles, &
                            shell_particles=shell_particles)
         CALL write_trajectory(force_env, root_section, it=it, time=0.0_dp, dtime=0.0_dp, &
                               etot=etot, pk_name="SHELL_TRAJECTORY", middle_name="shpos", &
                               particles=shell_particles)
         IF (shell_adiabatic) THEN
            CALL write_trajectory(force_env, root_section, it=it, time=0.0_dp, dtime=0.0_dp, &
                                  etot=etot, pk_name="SHELL_FORCES", middle_name="shfrc", &
                                  particles=shell_particles)
            CALL write_trajectory(force_env, root_section, it=it, time=0.0_dp, dtime=0.0_dp, &
                                  etot=etot, pk_name="CORE_TRAJECTORY", middle_name="copos", &
                                  particles=core_particles)
            CALL write_trajectory(force_env, root_section, it=it, time=0.0_dp, dtime=0.0_dp, &
                                  etot=etot, pk_name="CORE_FORCES", middle_name="cofrc", &
                                  particles=core_particles)
         END IF
      END IF

   END SUBROUTINE write_geo_traj

! **************************************************************************************************
!> \brief ...
!> \param gopt_env ...
!> \param output_unit ...
!> \param label ...
!> \date  01.2008
!> \author Teodoro Laino [tlaino] - University of Zurich
! **************************************************************************************************
   SUBROUTINE print_geo_opt_header(gopt_env, output_unit, label)

      TYPE(gopt_f_type), POINTER                         :: gopt_env
      INTEGER, INTENT(IN)                                :: output_unit
      CHARACTER(LEN=*), INTENT(IN)                       :: label

      CHARACTER(LEN=default_string_length)               :: my_format, my_label
      INTEGER                                            :: ix

      IF (output_unit > 0) THEN
         WRITE (UNIT=output_unit, FMT="(/,T2,A)") REPEAT("*", 79)
         IF (gopt_env%dimer_rotation) THEN
            my_label = "OPTIMIZING DIMER ROTATION"
         ELSE
            my_label = "STARTING "//gopt_env%tag(1:8)//" OPTIMIZATION"
         END IF

         ix = (80 - 7 - LEN_TRIM(my_label))/2
         ix = ix + 5
         my_format = "(T2,A,T"//cp_to_string(ix)//",A,T78,A)"
         WRITE (UNIT=output_unit, FMT=TRIM(my_format)) "***", TRIM(my_label), "***"

         ix = (80 - 7 - LEN_TRIM(label))/2
         ix = ix + 5
         my_format = "(T2,A,T"//cp_to_string(ix)//",A,T78,A)"
         WRITE (UNIT=output_unit, FMT=TRIM(my_format)) "***", TRIM(label), "***"

         WRITE (UNIT=output_unit, FMT="(T2,A)") REPEAT("*", 79)
         CALL m_flush(output_unit)
      END IF
   END SUBROUTINE print_geo_opt_header

! **************************************************************************************************
!> \brief ...
!> \param gopt_env ...
!> \param output_unit ...
!> \date  01.2008
!> \author Teodoro Laino [tlaino] - University of Zurich
! **************************************************************************************************
   SUBROUTINE print_geo_opt_nc(gopt_env, output_unit)

      TYPE(gopt_f_type), POINTER                         :: gopt_env
      INTEGER, INTENT(IN)                                :: output_unit

      IF (output_unit > 0) THEN
         WRITE (UNIT=output_unit, FMT="(/,T2,A)") &
            "*** MAXIMUM NUMBER OF OPTIMIZATION STEPS REACHED ***"
         IF (.NOT. gopt_env%dimer_rotation) THEN
            WRITE (UNIT=output_unit, FMT="(T2,A)") &
               "***        EXITING GEOMETRY OPTIMIZATION         ***"
         ELSE
            WRITE (UNIT=output_unit, FMT="(T2,A)") &
               "***        EXITING ROTATION OPTIMIZATION         ***"
         END IF
         CALL m_flush(output_unit)
      END IF

   END SUBROUTINE print_geo_opt_nc

! **************************************************************************************************
!> \brief   Prints information during GEO_OPT common to all optimizers
!> \param force_env ...
!> \param root_section ...
!> \param motion_section ...
!> \param its ...
!> \param opt_energy ...
!> \date    02.2008
!> \author  Teodoro Laino [tlaino] - University of Zurich
!> \version 1.0
! **************************************************************************************************
   SUBROUTINE geo_opt_io(force_env, root_section, motion_section, its, opt_energy)

      TYPE(force_env_type), POINTER                      :: force_env
      TYPE(section_vals_type), POINTER                   :: root_section, motion_section
      INTEGER, INTENT(IN)                                :: its
      REAL(KIND=dp), INTENT(IN)                          :: opt_energy

      TYPE(atomic_kind_list_type), POINTER               :: atomic_kinds
      TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
      TYPE(cell_type), POINTER                           :: cell
      TYPE(cp_subsys_type), POINTER                      :: subsys
      TYPE(distribution_1d_type), POINTER                :: local_particles
      TYPE(mp_para_env_type), POINTER                    :: para_env
      TYPE(particle_list_type), POINTER                  :: particles
      TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
      TYPE(virial_type), POINTER                         :: virial

      NULLIFY (para_env, atomic_kind_set, subsys, particle_set, &
               local_particles, atomic_kinds, particles)

      ! Write Restart File
      CALL write_restart(force_env=force_env, root_section=root_section)

      ! Write Trajectory
      CALL write_geo_traj(force_env, root_section, its, opt_energy)

      ! Write the stress Tensor
      CALL force_env_get(force_env, cell=cell, para_env=para_env, &
                         subsys=subsys)
      CALL cp_subsys_get(subsys=subsys, atomic_kinds=atomic_kinds, local_particles=local_particles, &
                         particles=particles, virial=virial)
      atomic_kind_set => atomic_kinds%els
      particle_set => particles%els
      CALL virial_evaluate(atomic_kind_set, particle_set, local_particles, &
                           virial, para_env)
      CALL write_stress_tensor(virial, cell, motion_section, its, 0.0_dp)

      ! Write the cell
      CALL write_simulation_cell(cell, motion_section, its, 0.0_dp)

   END SUBROUTINE geo_opt_io

! **************************************************************************************************
!> \brief   Apply coordinate transformations after cell (shape) change
!> \param gopt_env ...
!> \param cell ...
!> \param x ...
!> \param update_forces ...
!> \date    05.11.2012 (revised version of unbiase_coordinates moved here, MK)
!> \author  Matthias Krack
!> \version 1.0
! **************************************************************************************************
   SUBROUTINE apply_cell_change(gopt_env, cell, x, update_forces)

      TYPE(gopt_f_type), POINTER                         :: gopt_env
      TYPE(cell_type), POINTER                           :: cell
      REAL(KIND=dp), DIMENSION(:), POINTER               :: x
      LOGICAL, INTENT(IN)                                :: update_forces

      INTEGER                                            :: i, iatom, idg, j, natom, nparticle, &
                                                            shell_index
      REAL(KIND=dp)                                      :: fc, fs, mass
      REAL(KIND=dp), DIMENSION(3)                        :: s
      TYPE(cell_type), POINTER                           :: cell_ref
      TYPE(cp_subsys_type), POINTER                      :: subsys
      TYPE(particle_list_type), POINTER                  :: core_particles, particles, &
                                                            shell_particles

      NULLIFY (cell_ref)
      NULLIFY (core_particles)
      NULLIFY (particles)
      NULLIFY (shell_particles)
      NULLIFY (subsys)

      natom = force_env_get_natom(gopt_env%force_env)
      nparticle = force_env_get_nparticle(gopt_env%force_env)
      CALL force_env_get(gopt_env%force_env, &
                         subsys=subsys)
      CALL cp_subsys_get(subsys=subsys, &
                         core_particles=core_particles, &
                         particles=particles, &
                         shell_particles=shell_particles)

      ! Retrieve the reference cell
      CALL cell_create(cell_ref)
      CALL cell_copy(cell, cell_ref, tag="CELL_OPT_REF")

      ! Load the updated cell information
      SELECT CASE (gopt_env%cell_method_id)
      CASE (default_cell_direct_id)
         idg = 3*nparticle
         CALL init_cell(cell_ref, hmat=gopt_env%h_ref)
      CASE (default_cell_geo_opt_id, default_cell_md_id)
         idg = 0
      END SELECT
      CPASSERT((SIZE(x) == idg + 6))

      IF (update_forces) THEN

         ! Transform particle forces back to reference cell
         idg = 1
         DO iatom = 1, natom
            CALL real_to_scaled(s, x(idg:idg + 2), cell)
            CALL scaled_to_real(x(idg:idg + 2), s, cell_ref)
            idg = idg + 3
         END DO

      ELSE

         ! Update cell
         DO i = 1, 3
            DO j = 1, i
               idg = idg + 1
               cell%hmat(j, i) = x(idg)
            END DO
         END DO
         CALL init_cell(cell)
         CALL cp_subsys_set(subsys, cell=cell)

         ! Retrieve particle coordinates for the current cell
         SELECT CASE (gopt_env%cell_method_id)
         CASE (default_cell_direct_id)
            idg = 1
            DO iatom = 1, natom
               CALL real_to_scaled(s, x(idg:idg + 2), cell_ref)
               shell_index = particles%els(iatom)%shell_index
               IF (shell_index == 0) THEN
                  CALL scaled_to_real(particles%els(iatom)%r, s, cell)
               ELSE
                  CALL scaled_to_real(core_particles%els(shell_index)%r, s, cell)
                  i = 3*(natom + shell_index - 1) + 1
                  CALL real_to_scaled(s, x(i:i + 2), cell_ref)
                  CALL scaled_to_real(shell_particles%els(shell_index)%r, s, cell)
                  ! Update atomic position due to core and shell motion
                  mass = particles%els(iatom)%atomic_kind%mass
                  fc = core_particles%els(shell_index)%atomic_kind%shell%mass_core/mass
                  fs = shell_particles%els(shell_index)%atomic_kind%shell%mass_shell/mass
                  particles%els(iatom)%r(1:3) = fc*core_particles%els(shell_index)%r(1:3) + &
                                                fs*shell_particles%els(shell_index)%r(1:3)
               END IF
               idg = idg + 3
            END DO
         CASE (default_cell_geo_opt_id, default_cell_md_id)
            DO iatom = 1, natom
               shell_index = particles%els(iatom)%shell_index
               IF (shell_index == 0) THEN
                  CALL real_to_scaled(s, particles%els(iatom)%r, cell_ref)
                  CALL scaled_to_real(particles%els(iatom)%r, s, cell)
               ELSE
                  CALL real_to_scaled(s, core_particles%els(shell_index)%r, cell_ref)
                  CALL scaled_to_real(core_particles%els(shell_index)%r, s, cell)
                  i = 3*(natom + shell_index - 1) + 1
                  CALL real_to_scaled(s, shell_particles%els(shell_index)%r, cell_ref)
                  CALL scaled_to_real(shell_particles%els(shell_index)%r, s, cell)
                  ! Update atomic position due to core and shell motion
                  mass = particles%els(iatom)%atomic_kind%mass
                  fc = core_particles%els(shell_index)%atomic_kind%shell%mass_core/mass
                  fs = shell_particles%els(shell_index)%atomic_kind%shell%mass_shell/mass
                  particles%els(iatom)%r(1:3) = fc*core_particles%els(shell_index)%r(1:3) + &
                                                fs*shell_particles%els(shell_index)%r(1:3)
               END IF
            END DO
         END SELECT

      END IF

      CALL cell_release(cell_ref)

   END SUBROUTINE apply_cell_change

END MODULE gopt_f_methods
